/* Copyright (c) Colorado School of Mines, 2011.*/
/* All rights reserved.                       */

/*********************** self documentation **********************/
/*****************************************************************************
AIRY - Approximate the Airy functions  Ai(x), Bi(x) and their respective
       derivatives  Ai'(x), Bi'(x)

airya		return approximation of Ai(x)
airypa		return approximation of Ai'(x)
airyb		return approximation of Bi(x)
airybp		return approximation of Bi'(x)

******************************************************************************
Function Prototypes:
float airya (float x);
float airyap (float x);
float airyb (float x);
float airybp (float x);

******************************************************************************
Input:
x		value at which to evaluate Ai(x)

Returned:
airya		Ai(x)
airypa		Ai'(x)
airyb		Bi(x)
airybp		Bi'(x)

******************************************************************************
Reference:
The approximation is derived from tables and formulas in Abramowitz
and Stegun, p. 475-477.

******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/06/89
*****************************************************************************/
/**************** end self doc ********************************/
#include <cwp.h>

float airya (float x)
/*****************************************************************************
Return approximation to the Airy function Ai(x)
******************************************************************************
Input:
x		value at which to evaluate Ai(x)

Returned:	Ai(x)
******************************************************************************
Reference:
The approximation is derived from tables and formulas in Abramowitz
and Stegun, p. 475-477.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/06/89
*****************************************************************************/
{
	static float ap[101] = {
		0.35502, 0.35243, 0.34985, 0.34726, 0.34467,
		0.34209, 0.33951, 0.33693, 0.33435, 0.33177,
		0.32920, 0.32663, 0.32406, 0.32150, 0.31894,
		0.31639, 0.31384, 0.31130, 0.30876, 0.30623,
		0.30370, 0.30118, 0.29866, 0.29615, 0.29365,
		0.29116, 0.28867, 0.28619, 0.28372, 0.28126,
		0.27880, 0.27635, 0.27392, 0.27149, 0.26906,
		0.26665, 0.26425, 0.26186, 0.25947, 0.25710,
		0.25474, 0.25238, 0.25004, 0.24771, 0.24539,
		0.24308, 0.24078, 0.23849, 0.23621, 0.23394,
		0.23169, 0.22945, 0.22721, 0.22499, 0.22279,
		0.22059, 0.21841, 0.21624, 0.21408, 0.21193,
		0.20980, 0.20767, 0.20556, 0.20347, 0.20138,
		0.19931, 0.19726, 0.19521, 0.19318, 0.19116,
		0.18916, 0.18717, 0.18519, 0.18322, 0.18127,
		0.17933, 0.17741, 0.17549, 0.17360, 0.17171,
		0.16984, 0.16798, 0.16614, 0.16431, 0.16249,
		0.16069, 0.15890, 0.15713, 0.15536, 0.15362,
		0.15188, 0.15016, 0.14845, 0.14676, 0.14508,
		0.14341, 0.14176, 0.14012, 0.13850, 0.13689, 0.13529 };
	static float fp1[11] = {
		0.527027, 0.528783, 0.530601, 0.532488, 0.534448,
		0.536489, 0.538618, 0.540844, 0.543180, 0.545636, 0.548230 };
	static float fp2[11] = {
		0.548230, 0.549584, 0.550980, 0.552421, 0.553912,
		0.555456, 0.557058, 0.558724, 0.560462, 0.562280, 0.564190 };
	static float an1[101] = {
		0.35502, 0.35761, 0.36020, 0.36279, 0.36537,
		0.36796, 0.37054, 0.37312, 0.37560, 0.37827,
		0.38084, 0.38341, 0.38597, 0.38853, 0.39109,
		0.39364, 0.39618, 0.39871, 0.40124, 0.40376,
		0.40628, 0.40879, 0.41128, 0.41377, 0.41625,
		0.41872, 0.42118, 0.42363, 0.42606, 0.42849,
		0.43090, 0.43330, 0.43568, 0.43805, 0.44041,
		0.44275, 0.44508, 0.44739, 0.44968, 0.45196,
		0.45422, 0.45646, 0.45868, 0.46089, 0.46307,
		0.46523, 0.46738, 0.46950, 0.47159, 0.47367,
		0.47572, 0.47775, 0.47976, 0.48174, 0.48369,
		0.48562, 0.48752, 0.48939, 0.49124, 0.49306,
		0.49484, 0.49660, 0.49833, 0.50003, 0.50170,
		0.50333, 0.50493, 0.50650, 0.50803, 0.50953,
		0.51100, 0.51242, 0.51382, 0.51517, 0.51649,
		0.51777, 0.51901, 0.52021, 0.52137, 0.52249,
		0.52357, 0.52461, 0.52560, 0.52655, 0.52746,
		0.52832, 0.52914, 0.52991, 0.53064, 0.53132,
		0.53195, 0.53254, 0.53308, 0.53356, 0.53400,
		0.53439, 0.53473, 0.53501, 0.53525, 0.53543, 0.53556 };
	static float an2[91] = {
		 0.53556, 0.53381, 0.52619, 0.51227, 0.49170,
		 0.46425, 0.42986, 0.38860, 0.34076, 0.28680,
		 0.22740, 0.16348, 0.09614, 0.02670,-0.04333,
		-0.11232,-0.17850,-0.24003,-0.29509,-0.34190,
		-0.37881,-0.40438,-0.41744,-0.41718,-0.40319,
		-0.37553,-0.33477,-0.28201,-0.21885,-0.14741,
		-0.07026, 0.00967, 0.08921, 0.16499, 0.23370,
		 0.29215, 0.33749, 0.36736, 0.38003, 0.37453,
		 0.35076, 0.30952, 0.25258, 0.18256, 0.10293,
		 0.01778,-0.06833,-0.15062,-0.22435,-0.28512,
		-0.32914,-0.35351,-0.35652,-0.33734,-0.29713,
		-0.23802,-0.16352,-0.07831, 0.01210, 0.10168,
		 0.18428, 0.25403, 0.30585, 0.33577, 0.34132,
		 0.32177, 0.27825, 0.21372, 0.13285, 0.04170,
		-0.05270,-0.14290,-0.22159,-0.28223,-0.31959,
		-0.33029,-0.31311,-0.26920,-0.20205,-0.11726,
		-0.02213, 0.07495, 0.16526, 0.24047, 0.29347,
		 0.31910, 0.31465, 0.28023, 0.21886, 0.13623, 0.04024 };
	static float fn1[6] = {
		0.39752, 0.39781, 0.39809, 0.39838, 0.39866, 0.39894 };
	static float fn2[6] = {
		0.40028, 0.40002, 0.39975, 0.39949, 0.39921, 0.39894 };

	int ix,iy;
	float ax,frac,ai,sx,y,ay,f,f1,f2,oy;

	if (x>=0.0) {
		if (x<1.0) {
			ax = x*100.0;
			ix = ax;
			frac = ax-ix;
			ai = (1.0-frac)*ap[ix]+frac*ap[ix+1];
		} else {
			sx = sqrt(x);
			y = 1.5/(x*sx);
			if (y>0.5) {
				ay = 15.0-y*10.0;
				iy = ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				f = (1.0-frac)*fp1[iy]+frac*fp1[iy+1];
				ai = 0.5*exp(-1.0/y)*f/sqrt(sx);
			} else {
				ay = 10.0-y*20.0;
				iy = ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				f = (1.0-frac)*fp2[iy]+frac*fp2[iy+1];
				ai = 0.5*exp(-1.0/y)*f/sqrt(sx);
			}
		}
	} else {
		if (x>-10.0) {
			if (x>-1.0) {
				ax = -x*100.0;
				ix = ax;
				frac = ax-ix;
				ai = (1.0-frac)*an1[ix]+frac*an1[ix+1];
			} else {
				ax = (-x-1.0)*10.0;
				ix = ax;
				frac = ax-ix;
				ai = (1.0-frac)*an2[ix]+frac*an2[ix+1];
			}
		} else {
			sx = sqrt(-x);
			y = 1.5/(-x*sx);
			ay = 5.0-y*100.0;
			iy = ay;
			frac = ay-iy;
			f1 = (1.0-frac)*fn1[iy]+frac*fn1[iy+1];
			f2 = (1.0-frac)*fn2[iy]+frac*fn2[iy+1];
			oy = 1.0/y;
			ai = (f1*cos(oy)+f2*sin(oy))/sqrt(sx);
		}
	}
	return ai;
}


float airyap (float x)
/*****************************************************************************
 Return approximation to the derivative of the Airy function Ai'(x)
******************************************************************************
Input:
x		value at which to evaluate Ai'(x)

Returned:	Ai'(x)
******************************************************************************
Notes:
The approximation is derived from tables and formulas in Abramowitz
and Stegun, p. 475-477.
******************************************************************************
Authors:  Dave Hale, Craig Artley, Colorado School of Mines, 11/13/90
*****************************************************************************/
{
	static float app[101] = {
		-0.25881, -0.25880, -0.25874, -0.25866, -0.25854,
		-0.25838, -0.25819, -0.25797, -0.25772, -0.25744,
		-0.25713, -0.25678, -0.25641, -0.25600, -0.25557,
		-0.25511, -0.25462, -0.25411, -0.25356, -0.25300,
		-0.25240, -0.25178, -0.25114, -0.25047, -0.24977,
		-0.24906, -0.24832, -0.24756, -0.24677, -0.24597,
		-0.24514, -0.24429, -0.24343, -0.24254, -0.24164,
		-0.24071, -0.23977, -0.23881, -0.23783, -0.23684,
		-0.23583, -0.23480, -0.23376, -0.23270, -0.23163,
		-0.23054, -0.22944, -0.22833, -0.22720, -0.22606,
		-0.22491, -0.22374, -0.22257, -0.22138, -0.22018,
		-0.21897, -0.21775, -0.21653, -0.21529, -0.21404,
		-0.21279, -0.21153, -0.21025, -0.20898, -0.20769,
		-0.20640, -0.20510, -0.20380, -0.20248, -0.20117,
		-0.19985, -0.19852, -0.19719, -0.19585, -0.19451,
		-0.19317, -0.19182, -0.19047, -0.18912, -0.18777,
		-0.18641, -0.18505, -0.18369, -0.18232, -0.18096,
		-0.17959, -0.17823, -0.17686, -0.17549, -0.17413,
		-0.17276, -0.17139, -0.17003, -0.16866, -0.16730,
		-0.16593, -0.16457, -0.16321, -0.16185, -0.16050, -0.15914 };
	static float gp1[11] = {
		0.619954, 0.617156, 0.614275, 0.611305, 0.608239,
		0.605068, 0.601782, 0.598372, 0.594823, 0.591120, 0.587245 };
	static float gp2[11] = {
		0.587245, 0.585235, 0.583174, 0.581056, 0.578878,
		0.576635, 0.574320, 0.571927, 0.569448, 0.566873, 0.564190 };
	static float apn1[101] = {
		-0.25881, -0.25880, -0.25874, -0.25865, -0.25852,
		-0.25836, -0.25816, -0.25792, -0.25763, -0.25731,
		-0.25695, -0.25655, -0.25661, -0.25563, -0.25510,
		-0.25453, -0.25392, -0.25326, -0.25256, -0.25182,
		-0.25103, -0.25019, -0.24931, -0.24838, -0.24741,
		-0.24638, -0.24531, -0.24419, -0.24303, -0.24181,
		-0.24054, -0.23922, -0.23785, -0.23643, -0.23496,
		-0.23344, -0.23186, -0.23023, -0.22855, -0.22682,
		-0.22503, -0.22318, -0.22128, -0.21933, -0.21732,
		-0.21525, -0.21313, -0.21095, -0.20872, -0.20643,
		-0.20408, -0.20167, -0.19920, -0.19668, -0.19410,
		-0.19146, -0.18875, -0.18600, -0.18318, -0.18030,
		-0.17736, -0.17436, -0.17130, -0.16818, -0.16500,
		-0.16176, -0.15846, -0.15509, -0.15167, -0.14818,
		-0.14464, -0.14103, -0.13736, -0.13363, -0.12984,
		-0.12599, -0.12207, -0.11810, -0.11406, -0.10996,
		-0.10580, -0.10159, -0.09731, -0.09297, -0.08857,
		-0.08410, -0.07958, -0.07500, -0.07036, -0.06566,
		-0.06091, -0.05609, -0.05121, -0.04628, -0.04129,
		-0.03624, -0.03114, -0.02597, -0.02076, -0.01548, -0.01016 };
	static float apn2[91] = {
		-0.01016, 0.04602, 0.10703, 0.17199, 0.23981,
		 0.30918, 0.37854, 0.44612, 0.50999, 0.56809,
		 0.61825, 0.65834, 0.68624, 0.70003, 0.69801,
		 0.67885, 0.64163, 0.58600, 0.51221, 0.42118,
		 0.31458, 0.19482, 0.06503,-0.07096,-0.20874,
		-0.34344,-0.46986,-0.58272,-0.67688,-0.74755,
		-0.79062,-0.80287,-0.78221,-0.72794,-0.64085,
		-0.52336,-0.37953,-0.21499,-0.03676, 0.14695,
		 0.32719, 0.49458, 0.63990, 0.75457, 0.83122,
		 0.86419, 0.85003, 0.78781, 0.67943, 0.52962,
		 0.34593, 0.13836,-0.08106,-0.29899,-0.50147,
		-0.67495,-0.80711,-0.88790,-0.91030,-0.87103,
		-0.77100,-0.61552,-0.41412,-0.18009, 0.07027,
		 0.31880, 0.54671, 0.73605, 0.87115, 0.94004,
		 0.93556, 0.85621, 0.70659, 0.49727, 0.24422,
		-0.03231,-0.30933,-0.56297,-0.77061,-0.91289,
		-0.97566,-0.95149,-0.84067,-0.65149,-0.39986,
		-0.10809, 0.19695, 0.48628, 0.73154, 0.90781, 0.99626 };
	static float gn1[6] = {
		0.40092, 0.40052, 0.40012, 0.39972, 0.39933, 0.39894 };
	static float gn2[6] = {
		0.39704, 0.39741, 0.39779, 0.39817, 0.39855, 0.39894 };

	int ix,iy;
	float ax,frac,aip,sx,y,ay,g,g1,g2,oy;

	if (x>=0.0) {
		if (x<1.0) {
			ax = x*100.0;
			ix = ax;
			frac = ax-ix;
			aip = (1.0-frac)*app[ix]+frac*app[ix+1];
		} else {
			sx = sqrt(x);
			y = 1.5/(x*sx);
			if (y>0.5) {
				ay = 15.0-y*10.0;
				iy = ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				g = (1.0-frac)*gp1[iy]+frac*gp1[iy+1];
				aip = -0.5*exp(-1.0/y)*g*sqrt(sx);
			} else {
				ay = 10.0-y*20.0;
				iy = ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				g = (1.0-frac)*gp2[iy]+frac*gp2[iy+1];
				aip = -0.5*exp(-1.0/y)*g*sqrt(sx);
			}
		}
	} else {
		if (x>-10.0) {
			if (x>-1.0) {
				ax = -x*100.0;
				ix = ax;
				frac = ax-ix;
				aip = (1.0-frac)*apn1[ix]+frac*apn1[ix+1];
			} else {
				ax = (-x-1.0)*10.0;
				ix = ax;
				frac = ax-ix;
				aip = (1.0-frac)*apn2[ix]+frac*apn2[ix+1];
			}
		} else {
			sx = sqrt(-x);
			y = 1.5/(-x*sx);
			ay = 5.0-y*100.0;
			iy = ay;
			frac = ay-iy;
			g1 = (1.0-frac)*gn1[iy]+frac*gn1[iy+1];
			g2 = (1.0-frac)*gn2[iy]+frac*gn2[iy+1];
			oy = 1.0/y;
			aip = (g1*sin(oy)-g2*cos(oy))*sqrt(sx);
		}
	}
	return aip;
}


float airyb (float x)
/*****************************************************************************
Return approximation to the Airy function Bi(x)
******************************************************************************
Input:
x		value at which to evaluate Bi(x)

Returned:	Bi(x)
******************************************************************************
Notes:
The approximation is derived from tables and formulas in Abramowitz
and Stegun, p. 475-477.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 07/05/89
*****************************************************************************/
{
	static float bp[101] = {
		0.61492, 0.61940, 0.62389, 0.62837, 0.63286,
		0.63735, 0.64184, 0.64634, 0.65084, 0.65534,
		0.65986, 0.66438, 0.66890, 0.67343, 0.67798,
		0.68253, 0.68709, 0.69167, 0.69625, 0.70085,
		0.70546, 0.71008, 0.71472, 0.71938, 0.72405,
		0.72874, 0.73345, 0.73818, 0.74292, 0.74769,
		0.75248, 0.75729, 0.76213, 0.76699, 0.77187,
		0.77678, 0.78172, 0.78669, 0.79169, 0.79671,
		0.80177, 0.80686, 0.81198, 0.81714, 0.82233,
		0.82755, 0.83282, 0.83812, 0.84346, 0.84885,
		0.85427, 0.85974, 0.86525, 0.87081, 0.87641,
		0.88206, 0.88776, 0.89350, 0.89930, 0.90515,
		0.91106, 0.91702, 0.92303, 0.92910, 0.93524,
		0.94143, 0.94768, 0.95399, 0.96037, 0.96681,
		0.97332, 0.97990, 0.98655, 0.99327, 1.00006,
		1.00693, 1.01387, 1.02088, 1.02798, 1.03516,
		1.04242, 1.04976, 1.05719, 1.06470, 1.07230,
		1.07999, 1.08778, 1.09566, 1.10363, 1.11170,
		1.11987, 1.12814, 1.13651, 1.14499, 1.15357,
		1.16226, 1.17106, 1.17998, 1.18901, 1.19815, 1.20742 };
	static float fp1[11] = {
		0.619912, 0.620335, 0.620327, 0.619799, 0.618649,
		0.616764, 0.614022, 0.610309, 0.605543, 0.599723, 0.593015 };
	static float fp2[11] = {
		0.593015, 0.589451, 0.585855, 0.582330, 0.578985,
		0.575908, 0.573135, 0.570636, 0.568343, 0.566204, 0.564190 };
	static float bn1[101] = {
		0.61492, 0.61044, 0.60596, 0.60147, 0.59698,
		0.59249, 0.58800, 0.58351, 0.57901, 0.57450,
		0.56999, 0.56548, 0.56096, 0.55643, 0.55189,
		0.54735, 0.54280, 0.53824, 0.53367, 0.52909,
		0.52450, 0.51990, 0.51529, 0.51067, 0.50604,
		0.50139, 0.49674, 0.49207, 0.48738, 0.48268,
		0.47797, 0.47325, 0.46851, 0.46375, 0.45898,
		0.45419, 0.44939, 0.44457, 0.43974, 0.43488,
		0.43002, 0.42513, 0.42023, 0.41531, 0.41037,
		0.40541, 0.40043, 0.39544, 0.39043, 0.38540,
		0.38035, 0.37528, 0.37019, 0.36508, 0.35996,
		0.35481, 0.34965, 0.34446, 0.33936, 0.33403,
		0.32879, 0.32352, 0.31824, 0.31294, 0.30761,
		0.30227, 0.29691, 0.29153, 0.28612, 0.28070,
		0.27526, 0.26980, 0.26432, 0.25883, 0.25331,
		0.24777, 0.24222, 0.23665, 0.23106, 0.22545,
		0.21982, 0.21418, 0.20852, 0.20284, 0.19714,
		0.19143, 0.18570, 0.17996, 0.17420, 0.16842,
		0.16263, 0.15683, 0.15101, 0.14518, 0.13933,
		0.13347, 0.12760, 0.12171, 0.11582, 0.10991, 0.10399 };
	static float bn2[91] = {
		 0.10399, 0.04432,-0.01582,-0.07576,-0.13472,
		-0.19178,-0.24596,-0.29620,-0.34140,-0.38046,
		-0.41230,-0.43590,-0.45036,-0.45492,-0.44905,
		-0.43242,-0.40500,-0.36709,-0.31929,-0.26258,
		-0.19828,-0.12807,-0.05390, 0.02196, 0.09710,
		 0.16893, 0.23486, 0.29235, 0.33904, 0.37289,
		 0.39223, 0.39593, 0.38346, 0.35494, 0.31122,
		 0.25387, 0.18514, 0.10794, 0.02570,-0.05774,
		-0.13836,-0.21208,-0.27502,-0.32371,-0.35531,
		-0.36781,-0.36017,-0.33245,-0.28589,-0.22282,
		-0.14669,-0.06182, 0.02679, 0.11373, 0.19354,
		 0.26101, 0.31159, 0.34172, 0.34908, 0.33283,
		 0.29376, 0.23425, 0.15821, 0.07087,-0.02159,
		-0.11246,-0.19493,-0.26267,-0.31030,-0.33387,
		-0.33125,-0.30230,-0.24904,-0.17550,-0.08751,
		 0.00775, 0.10235, 0.18820, 0.25778, 0.30483,
		 0.32494, 0.31603, 0.27858, 0.21570, 0.13293,
		 0.03778,-0.06091,-0.15379,-0.23186,-0.28738,-0.31467 };
	static float fn1[6] = {
		0.39752, 0.39781, 0.39809, 0.39838, 0.39866, 0.39894 };
	static float fn2[6] = {
		0.40028, 0.40002, 0.39975, 0.39949, 0.39921, 0.39894 };

	int ix,iy;
	float ax,frac,bi,sx,y,ay,f,f1,f2,oy;

	if (x>=0.0) {
		if (x<1.0) {
			ax = x*100.0;
			ix = ax;
			frac = ax-ix;
			bi = (1.0-frac)*bp[ix]+frac*bp[ix+1];
		} else {
			sx = sqrt(x);
			y = 1.5/(x*sx);
			if (y>0.5) {
				ay = 15.0-y*10.0;
				iy = ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				f = (1.0-frac)*fp1[iy]+frac*fp1[iy+1];
				bi = exp(1.0/y)*f/sqrt(sx);
			} else {
				ay = 10.0-y*20.0;
				iy = ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				f = (1.0-frac)*fp2[iy]+frac*fp2[iy+1];
				bi = exp(1.0/y)*f/sqrt(sx);
			}
		}
	} else {
		if (x>-10.0) {
			if (x>-1.0) {
				ax = -x*100.0;
				ix = ax;
				frac = ax-ix;
				bi = (1.0-frac)*bn1[ix]+frac*bn1[ix+1];
			} else {
				ax = (-x-1.0)*10.0;
				ix = ax;
				frac = ax-ix;
				bi = (1.0-frac)*bn2[ix]+frac*bn2[ix+1];
			}
		} else {
			sx = sqrt(-x);
			y = 1.5/(-x*sx);
			ay = 5.0-y*100.0;
			iy = ay;
			frac = ay-iy;
			f1 = (1.0-frac)*fn1[iy]+frac*fn1[iy+1];
			f2 = (1.0-frac)*fn2[iy]+frac*fn2[iy+1];
			oy = 1.0/y;
			bi = (f2*cos(oy)-f1*sin(oy))/sqrt(sx);
		}
	}
	return bi;
}


float airybp (float x)
/*****************************************************************************
Return approximation to the derivative of the Airy function Bi'(x)
******************************************************************************
Input:
x		value at which to evaluate Bi'(x)

Returned:	Bi'(x)
******************************************************************************
Notes:
The approximation is derived from tables and formulas in Abramowitz
and Stegun, p. 475-477.
******************************************************************************
Authors:  Dave Hale, Craig Artley, Colorado School of Mines, 11/13/90
*****************************************************************************/
{
	static float bpp[101] = {
		0.44828, 0.44831, 0.44841, 0.44856, 0.44878,
		0.44907, 0.44942, 0.44984, 0.45033, 0.45088,
		0.45151, 0.45220, 0.45297, 0.45381, 0.45472,
		0.45571, 0.45677, 0.45791, 0.45912, 0.46041,
		0.46178, 0.46324, 0.46477, 0.46638, 0.46808,
		0.46986, 0.47172, 0.47367, 0.47571, 0.47783,
		0.48004, 0.48235, 0.48474, 0.48722, 0.48980,
		0.49247, 0.49524, 0.49810, 0.50106, 0.50412,
		0.50728, 0.51053, 0.51389, 0.51736, 0.52092,
		0.52459, 0.52837, 0.53225, 0.53625, 0.54035,
		0.54457, 0.54890, 0.55334, 0.55789, 0.56257,
		0.56736, 0.57227, 0.57730, 0.58246, 0.58774,
		0.59314, 0.59867, 0.60433, 0.61012, 0.61603,
		0.62209, 0.62827, 0.63460, 0.64106, 0.64766,
		0.65440, 0.66129, 0.66832, 0.67549, 0.68282,
		0.69029, 0.69792, 0.70571, 0.71365, 0.72174,
		0.73000, 0.73842, 0.74701, 0.75576, 0.76468,
		0.77378, 0.78304, 0.79249, 0.80211, 0.81191,
		0.82190, 0.83207, 0.84243, 0.85298, 0.86373,
		0.87467, 0.88581, 0.89716, 0.90871, 0.92046, 0.93243 };
	static float gp1[11] = {
		0.478728, 0.479925, 0.481658, 0.484018, 0.487107,
		0.491037, 0.495921, 0.501859, 0.508909, 0.517032, 0.526011 };
	static float gp2[11] = {
		0.526011, 0.530678, 0.535345, 0.539902, 0.544235,
		0.548255, 0.551930, 0.555296, 0.558428, 0.561382, 0.564190 };
	static float bpn1[101] = {
		0.44828, 0.44831, 0.44841, 0.44856, 0.44877,
		0.44903, 0.44936, 0.44974, 0.45017, 0.45066,
		0.45121, 0.45180, 0.45245, 0.45315, 0.45390,
		0.45470, 0.45554, 0.45643, 0.45737, 0.45835,
		0.45938, 0.46045, 0.46156, 0.46272, 0.46391,
		0.46515, 0.46642, 0.46773, 0.46908, 0.47046,
		0.47188, 0.47333, 0.47481, 0.47632, 0.47787,
		0.47944, 0.48105, 0.48268, 0.48434, 0.48602,
		0.48773, 0.48946, 0.49122, 0.49299, 0.49479,
		0.49660, 0.49844, 0.50029, 0.50215, 0.50403,
		0.50593, 0.50784, 0.50976, 0.51169, 0.51363,
		0.51557, 0.51753, 0.51949, 0.52145, 0.52342,
		0.52540, 0.52737, 0.52934, 0.53132, 0.53329,
		0.53525, 0.53721, 0.53917, 0.54112, 0.54306,
		0.54499, 0.54692, 0.54883, 0.55072, 0.55260,
		0.55447, 0.55632, 0.55815, 0.55996, 0.56176,
		0.56353, 0.56527, 0.56699, 0.56869, 0.57036,
		0.57200, 0.57362, 0.57520, 0.57675, 0.57826,
		0.57974, 0.58119, 0.58260, 0.58397, 0.58530,
		0.58659, 0.58783, 0.58904, 0.59019, 0.59131, 0.59237 };
	static float bpn2[91] = {
		 0.59237, 0.60011, 0.60171, 0.59592, 0.58165,
		 0.55790, 0.52389, 0.47906, 0.42315, 0.35624,
		 0.27879, 0.19168, 0.09622,-0.00581,-0.11223,
		-0.22042,-0.32739,-0.42989,-0.52445,-0.60751,
		-0.67561,-0.72544,-0.75412,-0.75926,-0.73920,
		-0.69311,-0.62117,-0.52461,-0.40581,-0.26829,
		-0.11667, 0.04347, 0.20575, 0.36320, 0.50858,
		 0.63474, 0.73494, 0.80328, 0.83508, 0.82721,
		 0.77841, 0.68948, 0.56345, 0.40555, 0.22307,
		 0.02511,-0.17783,-0.37440,-0.55300,-0.70247,
		-0.81289,-0.87622,-0.88697,-0.84276,-0.74461,
		-0.59717,-0.40856,-0.19009, 0.04437, 0.27926,
		 0.49824, 0.68542, 0.82650, 0.90998, 0.92812,
		 0.87780, 0.76095, 0.58474, 0.36122, 0.10670,
		-0.15945,-0.41615,-0.64232,-0.81860,-0.92910,
		-0.96296,-0.91547,-0.78882,-0.59221,-0.34136,
		-0.05740, 0.23484, 0.50894, 0.73928, 0.90348,
		 0.98471, 0.97349, 0.86898, 0.67936, 0.42147, 0.11941 };
	static float gn1[6] = {
		0.40092, 0.40052, 0.40012, 0.39972, 0.39933, 0.39894 };
	static float gn2[6] = {
		0.39704, 0.39741, 0.39779, 0.39817, 0.39855, 0.39894 };

	int ix,iy;
	float ax,frac,bip,sx,y,ay,g,g1,g2,oy;

	if (x>=0.0) {
		if (x<1.0) {
			ax = x*100.0;
			ix = ax;
			frac = ax-ix;
			bip = (1.0-frac)*bpp[ix]+frac*bpp[ix+1];
		} else {
			sx = sqrt(x);
			y = 1.5/(x*sx);
			if (y>0.5) {
				ay = 15.0-y*10.0;
				iy = ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				g = (1.0-frac)*gp1[iy]+frac*gp1[iy+1];
				bip = exp(1.0/y)*g*sqrt(sx);
			} else {
				ay = 10.0-y*20.0;
				iy = ay;
				if (iy<0) iy = 0;
				else if (iy>9) iy = 9;
				frac = ay-iy;
				g = (1.0-frac)*gp2[iy]+frac*gp2[iy+1];
				bip = exp(1.0/y)*g*sqrt(sx);
			}
		}
	} else {
		if (x>-10.0) {
			if (x>-1.0) {
				ax = -x*100.0;
				ix = ax;
				frac = ax-ix;
				bip = (1.0-frac)*bpn1[ix]+frac*bpn1[ix+1];
			} else {
				ax = (-x-1.0)*10.0;
				ix = ax;
				frac = ax-ix;
				bip = (1.0-frac)*bpn2[ix]+frac*bpn2[ix+1];
			}
		} else {
			sx = sqrt(-x);
			y = 1.5/(-x*sx);
			ay = 5.0-y*100.0;
			iy = ay;
			frac = ay-iy;
			g1 = (1.0-frac)*gn1[iy]+frac*gn1[iy+1];
			g2 = (1.0-frac)*gn2[iy]+frac*gn2[iy+1];
			oy = 1.0/y;
			bip = (g1*cos(oy)+g2*sin(oy))*sqrt(sx);
		}
	}
	return bip;
}

#ifdef TEST
main()
{
	int nx,i;
	float xmin,xmax,x,dx;

	scanf("%d %*[^\n]", &nx);
	scanf("%f %*[^\n]", &xmin);
	scanf("%f %*[^\n]", &xmax);

	for (i=0,x=xmin,dx=(xmax-xmin)/(nx-1); i<nx; i++,x+=dx)
		printf("ai[%f] = %f\n", x, airya(x));
		printf("aip[%f] = %f\n", x, airyap(x));
		printf("bi[%f] = %f\n", x, airyb(x));
		printf("bip[%f] = %f\n", x, airybp(x));
}
#endif

